#####################################################################
# Results reverse causality 
#####################################################################

rm(list=ls())

library(lmtest)
library(sandwich)
library(stargazer)

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<40){
    warning("Fewer than 40 clusters, variances may be unreliable")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##############################################
# Read data and prepare basic output
##############################################

# Generate log  
#sink("008_outcome_rev_cau.txt")

# Read data 
load("006_rev_cau.Rdata")

##########################
# Regressions results
##########################

# Egotropic only treatment
mod1a = lm(outcome_ego ~ t_ind, data=d_match)
mod1 <- coeftest(mod1a, vcov = vcovCluster(mod1a, cluster = d_match$bairroibge))

# Egotropic controls and fixed effects
mod2a = lm(outcome_ego ~ t_ind + 
                         v5a_econpastego_imp + 
                         v5d_econpastbrazil_imp + 
                         v41e_thermserra_imp +
                         v50_ideology_imp + 
                         as.factor(bairroibge),
                         data = d_match)
mod2 <- coeftest(mod2a, vcov = vcovCluster(mod2a, cluster = d_match$bairroibge))

# Sociotropic only treatment
mod3a = lm(outcome_socio ~ t_ind, data=d_match)
mod3 <- coeftest(mod3a, vcov = vcovCluster(mod3a, cluster = d_match$bairroibge))

# Sociotropic controls and fixed effects
mod4a = lm(outcome_socio ~ t_ind + 
                           v5a_econpastego_imp + 
                           v5d_econpastbrazil_imp + 
                           v41e_thermserra_imp +
                           v50_ideology_imp +
                           as.factor(bairroibge),
                           data = d_match)
mod4 <- coeftest(mod4a, vcov = vcovCluster(mod2a, cluster = d_match$bairroibge))

# TABLE 3 MANUSCRIPT 
dim(d_match)
stargazer(mod2,mod4,
          type = "text",
          title="Regression results defection from the incumbent",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          #single.row = TRUE,
          #column.separate = c(4,4),
          column.labels   = c("Egotropic", "Sociotropic"),
          column.separate = c(1, 1),
          covariate.labels=c("Defection from the incumbent"),
          dep.var.caption  = "Economic evaluations",
          dep.var.labels.include = FALSE,
          omit = c("v5a_econpastego_imp", "v5d_econpastbrazil_imp","v41e_thermserra_imp","v50_ideology_imp","bairroibge","Constant"),
          add.lines = list(c("Controls","Yes","Yes"),
                           c("Neighborhood Fixed effects","Yes","Yes"),
                           c("Observations","330","330")))
          
#sink()
